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Neural firing is often subject to negative feedback by adaptation currents. These currents 
can induce strong correlations among the time intervals between spikes. Here we 
study analytically the interval correlations of a broad class of noisy neural oscillators 
with spike-triggered adaptation of arbitrary strength and time scale. Our weak-noise 
theory provides a general relation between the correlations and the phase-response 
curve (PRC) of the oscillator, proves anti-correlations between neighboring intervals for 
adapting neurons with type I PRC and identifies a single order parameter that determines 
the qualitative pattern of correlations. Monotonically decaying or oscillating correlation 
structures can be related to qualitatively different voltage traces after spiking, which can 
be explained by the phase plane geometry. At high firing rates, the long-term variability 
of the spike train associated with the cumulative interval correlations becomes small, 
independent of model details. Our results are verified by comparison with stochastic 
simulations of the exponential, leaky, and generalized integrate-and-fire models with 
adaptation. 



Keywords: spike-frequency adaptation, non-renewal process, serial correlation coefficient, phase-response curve, 
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1. INTRODUCTION 

The nerve cells of the brain are complex physical systems. They 
generate action potentials (spikes) by a nonlinear, adaptive, and 
noisy mechanism. In order to understand signal processing in 
single neurons, it is vital to analyze the sequence of the inter- 
spike intervals (ISIs) between adjacent action potentials. There 
is experimental evidence accumulating that the spiking in many 
cases is not a renewal process, i.e., a spike train with mutu- 
ally independent ISIs, but that intervals are typically correlated 
over a few lags (Lowen and Teich, 1992; Ratnam and Nelson, 
2000; Neiman and Russell, 2001; Nawrot et al., 2007; Engel et al., 
2008) [further reports are reviewed in (Farkhooi et al., 2009; 
Avila-Akerberg and Chacron, 2011)]. These correlations are a 
basic statistics of any spike train with important implications for 
information transmission and signal detection in neural systems 
(Ratnam and Nelson, 2000; Chacron et al, 2001, 2004; Avila- 
Akerberg and Chacron, 2011) and man-made signal detectors 
(Nikitin et al., 2012). They are often characterized by the serial 
correlation coefficient (SCC) 

( = {(T,-{m(T l + k -(T l + k ))) 
Pk ((T; - <r,»2> 

where T,- and T,' + j- are two ISIs lagged by an integer k and (•) 
denotes ensemble averaging. ISI correlations can be induced via 
correlated input to the neural dynamics, e.g. in the form of exter- 
nal colored noise (Middleton et al, 2003; Lindner, 2004), intrinsic 
noise from ion channels with slow kinetics (Fisch et al., 2012), or 
stochastic narrow-band input (Neiman and Russell, 2001, 2005; 
Bauermeister et al., 2013). 



Another ubiquitous mechanism for ISI correlations are 
slow feedback processes mediating spike-frequency adaptation 
(Chacron et al, 2000; Liu and Wang, 2001; Benda et al, 2005)— 
a phenomenon describing the reduced neuronal response to 
slowly changing stimuli (Benda and Herz, 2003; Gabbiani and 
Krapp, 2006). In the stationary state, these adaptation mech- 
anisms are typically associated with short-range correlations 
with a negative SCC at lag k = 1 and a reduced Fano factor 
as demonstrated by several numerical (Geisler and Goldberg, 
1966; Wang, 1998; Liu and Wang, 2001; Benda et al, 2010) 
and analytical studies (Schwalger et al., 2010; Schwalger and 
Lindner, 2010; Farkhooi et al, 2011; Urdapilleta, 2011). The cor- 
relation structure of adapting neurons can show qualitatively 
different patterns, ranging from monotonically decaying corre- 
lations to damped oscillations when plotted as a function of the 
lag (Ratnam and Nelson, 2000). Because ISI correlations shape 
spectral measures (Chacron et al., 2004), they bear implications 
for neural computation in general. However, a simple theory 
that predicts and explains possible correlation patterns is still 
lacking. 

In this article, we present a relation between the ISI correla- 
tion coefficient pj- and a basic characteristics of nonlinear neural 
dynamics, the phase-response curve (PRC). The PRC quantifies the 
advance (or delay) of the next spike caused by a small depolarizing 
current applied at a certain time after the last spike (Ermentrout, 
1996). For neurons which integrate up their input (integrator 
neurons), the PRC is positive at all times (type I PRC) whereas 
neurons, which show subthreshold resonances (resonator neu- 
rons), possess a PRC that is partly negative (type II PRC) 
(Ermentrout, 1996; Izhikevich, 2005; Ermentrout and Terman, 
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FIGURE 1 | Correlation patterns in the adaptive exponential IF model 
with x a = 10, y = 1, A 7 = 0.1, uj = 2 , D = 0.1. Adaptation is weak 
(A = 1, |i = 15) in (A) and strong (A = 10, \i = 80) in (B). Membrane 
voltage v(t) and adaptation variable a(f) with ISI sequences {T;J and peak 
adaptation values {a,) are shown in (A1.B1); time is in units of the 
membrane time constant x m . Colored pieces of trajectories in the phase 
plane (v, a) in (A2.B2) correspond to the respective colors in (A1.B1). The 
deterministic limit cycle (LC), determined by the initial (post-spike) values 
v = 0, a = a*, is indicated by a thick black line. For weak adaptation (A2) a 
short ISI Tj causes positive deviations 8a,- = a,- — a* and Sa, +1 = a, +1 — a* of 
peak values leading to long ISIs 7, +1 and 7/ + 2 and, hence, to a negative ISI 
correlation at all lags (A3). Because of the qualitatively different limit cycle 
for strong adaptation (B2), deviations 8a/ and 8a, +1 differ in sign, yielding an 
oscillatory correlation pattern (B3). 



2010). Below we show that resonator neurons possess a richer 
repertoire of correlation patterns than integrator neurons do. 

2. RESULTS 
2.1. MODEL 

Spike frequency adaptation can be modeled by Hodgkin-Huxley 
type neurons with a depolarization-activated adaptation current 
(Wang, 1998; Ermentrout et al., 2001; Benda and Herz, 2003). 
However, the spiking of such conductance-based models can in 
many instances be approximated by simpler multi-dimensional 
integrate-fire (IF) models that are equipped with a spike-triggered 
adaptation current (Treves, 1993; Izhikevich, 2003; Brette and 
Gerstner, 2005); adapting IF models perform excellently in 
predicting spike times of real cells under noisy stimulation 
(Gerstner and Naud, 2009). Here, we consider a stochastic non- 
linear multi-dimensional IF model for the membrane potential 
v, N auxiliary variables w; (j = 1, . . . , N) and a spike-triggered 
adaption current a(t): 

v=/o(v,w) + |A-a + i;(0, (2a) 
Wj=fj(v,vr), (2b) 

x a 'a = —a + T a A 8 (f - f,) . (2c) 
i 

The membrane potential v(t) is subject to weak Gaussian noise 
%(t) with (£(t)£(t0) = 2D8(f- t') and noise intensity D. The 
dynamics is complemented by a spike-and-reset mechanism: 
whenever V(t) reaches a threshold v(t), a spike is registered at 
time ti = t and vit) and w(f) = [wi(f), . . . , WN(t)] T are reset 
to v(fc^) = 0 and w(f + ) = w r (where t + denotes the right-sided 
limit t — > t j + 0). At the same time, a(t) suffers a jump by A > 
0 as seen from Equation (2c), which resembles high-threshold 
adaptation currents (Wang, 1998; Liu and Wang, 2001). The 
constant input current |x is assumed to be sufficiently large to 
ensure ongoing spiking even in the absence of noise. Note that 
the model is non-dimensionalized by measuring time in units of 
the membrane time constant x m ~ 10 ms and voltage in units of 
the distance between reset and spike-initiating potential (a typ- 
ical value is 15 mV). In particular, the adaptation time constant 
x a is measured relative to x m and the unit of the firing rate is 
t^ 1 ~ 100 Hz. 

An important special case, the adaptive exponential integrate- 
and-fire model (Brette and Gerstner, 2005) with purely spike- 
triggered adaptation and a white noise current with constant 
mean is illustrated in Figure 1. It assumes an exponential nonlin- 
earity /o(v) = — yv + yAj exp[(v — 1)/Aj] (Fourcaud-Trocme 
et al, 2003; Badel et al., 2008) and corresponds to N = 0. Time 
courses of v(t) and a(t) are shown in Figures 1A1.B1 for two dis- 
tinct correlation patterns possible in this model. The ISIs T; = 
ti — tj-\ are obtained as differences between subsequent spiking 
times f;. The sequence T;, T >+ i, T, + 2 displays patterns of short- 
long-long (Figure 1A1) and short-long-short (Figure 1B1), corre- 
sponding to a negative SCC, which decays monotonically with the 
lag k (Figure 1A3) or to an SCC oscillating with k (Figure 1B3). 
In the following, we develop a theory to analyze these and other 



correlation patterns possible in multi-dimensional adapting IF 
models. 

2.2. GENERAL THEORY 

In our model Equation (2), a(t) is the only variable that keeps a 
memory of the previous spike times thereby inducing correlations 
between ISIs. Over one ISI the time course of adaptation is an 
exponential decay, relating two adjacent peak values a, = 
and a; +1 = a{tf +l ) by 

a i+ i =a,e- T ' + l/t « + A (3) 

(Figures A 1,B1). We assume that in the deterministic case 
(D = 0) our model has a finite period T* (i.e., the model oper- 
ates in the tonically firing regime) and, hence, for D = 0 the map 
(3) has a stable fixed point 

a* = A/[l-exp(-r/t a )]. (4) 

The asymptotic deterministic dynamics can be interpreted as 
a limit-cycle like motion in the phase space from the reset 
point to the threshold and back by the instantaneous reset [cf. 
Figures 1A2,B2]. 

Weak noise will cause small deviations in the period 8T; = 
Ti — T* ~ Ti — (Tj) that are mutually correlated with coefficient 
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p k = (8T;8X,-_|_jt)/(8X 2 ). The peak adaptation values, however, 
also fluctuate, 8a; = a, — a*, and both deviations are related by 
linearizing Equation (3): 



8X, 



i+l 



2* V 



e r/Tfl 8a 



) 



(5) 



A second relation between the small deviations can be gained 
by considering how a small perturbation in the voltage dynam- 
ics affects the length of the period. This effect is captured by 
the infinitesimal phase response curve (PRC), Z(t), t e (0, T*) 
(Izhikevich, 2005; Ermentrout and Terman, 2010) (see Section 
4 for the precise definition). During the interval X;+i, the volt- 
age dynamics in Equation (2a) can be written as v = fo(v, w) + 
u, — (a* + hai)e~ ( - t ~ ti) t Xa + %(t). Compared to the deterministic 
limit cycle, the dynamics is perturbed by the weak noise and the 
small deviation in the adaptation Sflj-e - ^ - t "' x " yielding in linear 
response 

SX, + 1 = J dtZ(t) (hcne-i - t-(ti + 0) ■ (6) 
Combining Equations (5), (6) we obtain the stochastic map 



i+l = a#8a, + S,-, 



(7) 



where S; = ^r- f n °° df Z(t)%(t, + t) are uncorrected Gaussian 
random numbers and 



a* f 1 
1 / dtZ(t)e * . 

to JO 



(8) 



Note that local stability of the fixed point a* requires that |on>| < 
1. The covariance Cjfc = (8a;8a,-|-it) of the auto-regressive process 
Equation (7) can be calculated by elementary means and using 
Equation (5) we obtain for k > 1: 



Pit 



-A(l - ft) (a&) 



k- 1 



a(l - a 2 ft) 
1 + a 2 - 2a 2 ft' 



(9) 



In order to compute a and ft via Equation (8), we have to calculate 
X* and Z(t) (a* then follows from Equation (4)), which can be 
done for simple systems analytically. 

Our main result, Equations (8), (9), allows to draw a num- 
ber of general conclusions. It shows that the SCC is always a 
geometric sequence with respect to the lag k that can generate 
qualitatively different correlation patterns depending on the value 
of i"> and thus on PRC and adaptation current. Because |ai}| < 1 
and 0 < a < 1, the prefactor A in Equation (9) is always positive. 
Consequently, pi is negative for ft < 1 and positive ior ft > 1. 
Looking at Equation (8), we find that a positive PRC inevitably 
yields ft < 1 . This implies that adapting neurons with type I PRC 
possess negative correlations between adjacent ISIs. Intuitively, a 
short ISI causes in the following on average a higher inhibitory 
adaptation during the subsequent ISI. Such an inhibitory current 
always enlarges the ISI in type I neurons — hence, a short ISI is 
followed by a long ISI. 



The sign of higher lags is determined by the base of the power: 
for ft > 0 correlations decay monotonically, whereas for ft < 0 
the SCC oscillates. Two special cases are ft = 0 with a negative 
correlation at lag 1 and vanishing correlations at all higher lags 
and ft = 1 where all correlations vanish. Overall, we find five 
basic patterns corresponding to the cases — a -1 < ft < 0, ft = 0, 
0<i}<l,i}=l and 1 < ft < a -1 . These basic patterns cover 
all interval correlations discussed in previous theoretical studies 
(Schwalger and Lindner, 2010; Urdapilleta, 2011). Our geomet- 
ric formula generalizes the theory for the perfect IF model with 
adaptation (Schwalger et al., 2010) to more realistic, nonlinear 
multi-dimensional IF models with adaptation. 

The cumulative effect of the correlations can be described by 
the sum over all p^, which determines the long-time limit of the 
Fano factor and the low-frequency limit of the spike train power 
spectrum (for a definition of these quantities, see Section 4.2). 
Evaluating the geometric series yields 



k= 1 



A (1 - ft) 
1 - aft 



(10) 



This shows that adaptation in neurons with type I resetting 
(ft < 1) leads to a negative summed correlation and hence a 
reduced long-term variability. Furthermore, at high firing rates 
achieved by a strong input current u,, the sum in Equation (10) 
can be approximated by 



k=l 



+ 



1/2 



2 (1 + Ax a /v T ) z 



T* « X a . 



(ID 



In particular, for strong adaptation (Ax a 3> Vj) the sum is only 
slightly larger than —1/2. Note that by virtue of the funda- 
mental relation lim^oo P(f) = Cy (l + 2 J^T= l Pk) (Cox and 
Lewis, 1966) (see Section 4.2), the smallest possible value for the 
sum is — 1/2 in order to ensure the non-negativity of the Fano 
factor F(t). At this minimal value the long-term variability as 
expressed by the Fano factor vanishes even for a non-vanishing 
ISI variability as quantified by the coefficient of variation Cy. The 
latter quantity can also be estimated using the weak-noise the- 
ory: From Equation (7) one can calculate the variance of a; and 
using Equation (5) an approximation for Cy ~ (&T 2 )/T* 2 can be 
obtained as follows: 



2D- 



1 + a 2 - 2a l ft 



[1 



(aft) 2 ]T* 2 



L 



V 



df[Z(f)] 2 . 



(12) 



2.3. ONE-DIMENSIONAL IF MODELS WITH ADAPTATION 

In the simplest case (N = 0,/o(v, w) = f(v)) the PRC reads 

Z(f) = Z(T*) exp [J d//'(v 0 (f'))], 



(13) 



where Vn(f) is the limit cycle solution and Z(T*) = [f(vj) + u, — 
a* + A] -1 is the inverse of the velocity vn(X*) at the thresh- 
old, which is always positive. Thus, the PRC is positive for all 
t € (0, T*), i.e., one-dimensional IF models show type I behav- 
ior. From our general considerations, this implies a negative SCC 
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FIGURE 2 | ISI correlations and coefficient of variation (CV) of the 
adapting LIF model vs. firing rate 1/(7/) « 1/7*, where the rate is 
varied by increasing yi. The gray-shaded area corresponds to the 
fluctuation-driven regime (u, < yVf), where the assumptions of the theory 
do not hold. The panels display (from top to bottom) pi , p2, the sum 
5^™ =1 Pk and the CV for simulation (circles, m = 100) and theory (solid 
lines, m— > oo). (A) Moderate adaptation: A = 1, (B) strong adaptation: 
A = 10. Both: y = 1 , x a = 10, 0 = 0.1, Vr = 1 ■ Note that the firing rate is 
given in units of the inverse membrane time constant x^ 1 . 



at lag k = 1. The sign of the correlations at higher lags can be 
inferred from the sign of i}, for which one can show (Section 4) 
that 

1} = (f(0) + ii - a*)Z(0). (14) 

Because Z(0) > 0, the sign of t) is determined by the sign of 
/(0) + |x — a*. For weak adaptation such that a* </(0) + u, 
(achieved by a sufficiently small value of A or x a , Figure 1A), we 
will have $ > 0 and a negative correlation at all lags (Figure 1A3). 
In this case, a short ISI occurring by fluctuation will cause a pos- 
itive deviation 8a,- (Figure 1A2, green arrow). Geometrically, it is 
plausible that such a positive deviation causes a likewise positive 
deviation 8a, + i in the subsequent cycle (Figure 1A2, red arrow). 
Because a positive deviation is associated with a long ISI, the 
initial short ISI is on average followed by longer ISIs. 

In marked contrast, for strong adaptation such that a* > 
/(0) + u, (achieved by a sufficiently large value of A or t„), 
becomes negative and hence the SCC's sign alternates with the lag. 
This alternation of the sign can be understood by means of the 
phase plane. Let us again consider a positive deviation 8a, due to 
a short preceding ISI (Figure 1B2, green arrow). Because vn(O) = 
/(0) + [L — a* < 0, the neuron is reset above the v-nullcline and 
hence hyperpolarizes at the beginning of the interval, i.e., the 
trajectory makes a detour into the region of negative voltage (cor- 
responding to a "broad reset" in Naud et al. (2008)). A positive 
deviation 8a, leads to a larger detour (green trajectory) causing a 
sign inversion and hence a negative deviation 8a, +1 (Figure 1B2, 
red arrow). Because a positive (negative) deviation corresponds 
on average to a long (short) ISI, the alternation of 8a, also 
entails an alternation of the ISI correlations. Thus, the distinction 
between monotonic and alternating patterns relates to a qualita- 
tive distinction of the voltage trace after resetting [cf. "sharp" vs. 
"broad" resets in Naud et al. (2008)]. 

As demonstrated in Figures 1A3,B3, our theory works well 
for the adapting exponential integrate-and-fire model. We next 
demonstrate the validity of our approach over a broad range 
of firing rates (Figure 2) for another important ID model, the 
adapting leaky integrate-and-fire model (Treves, 1993) for which 
/(v) = — yv and 

Z(f) = exp[y(f - T*)]/(u. — yvj — a* + A) (15) 

(here T* has still to be determined from a transcendental 
equation). Changing the firing rate by varying the input cur- 
rent u,, we find a good agreement for the first two correlation 
coefficients and the sum of all p^; the approximation of the CV 
shows deviations from simulation results when the input current 
[L becomes small (approaching the fluctuation-driven regime). In 
accordance with previous findings (Wang, 1998; Liu and Wang, 
2001; Benda et al, 2010; Nesse et al, 2010; Schwalger et al, 
2010; Schwalger and Lindner, 2010; Urdapilleta, 2011), the first 
correlation coefficient pi displays a minimum corresponding to 
strong anti-correlations between adjacent intervals. The correla- 
tions at lag 2 can be positive for a finite range of firing rates if 
the adaptation strength is sufficiently large (Figure 2B), whereas 
for moderate adaptation we find a negative p2 at all firing rates 



(Figure 2A). In both cases, however, the sum of SCCs approaches 
a value close to — 1 /2 for high firing rates as predicted by Equation 
(11) (Figure 2, bottom). This is strikingly similar to experimen- 
tal data from weakly electric fish, in which some electro-receptors 
display a monotonically decaying SCC and some show an oscil- 
latory SCC (Ratnam and Nelson, 2000) but all cells exhibit a 
sum close to —1/2 (Ratnam and Goense, 2004). Finally, Figure 2 
reveals a local maximum of the CV for some suprathreshold 
current u, — an effect that has been described by Nesse et al. 
(2008). 

2.4. GENERALIZED INTEGRATE-AND-FIRE MODEL WITH ADAPTATION 

Different correlation patterns become possible if we consider a 
type II PRC, which is by definition partly negative and can lead 
to a negative value of the integral in Equation (8), and hence to 
i} > 1. This corresponds to a non-negative SCC at lag 1, which 
is infeasible in the one-dimensional case. To test the prediction 
Pi > 0, we study the generalized integrate-and-fire (GIF) model 
(Brunei et al, 2003) with spike-triggered adaptation. This model 
is defined by /o(v, w) = —yv — pV and /i(v, w) = (v — w)/x w . 
Using the method described in Section 4, the PRC is obtained as 

ez (f - r) Tcos(f2(f - T*)) - sin (C2(f - T*))l 

Z(f) = ^ i (16) 

u, — yvj — pVoCT*) — a* + A 



where v = y + l/x w , £2 = J — ~ and Wo(t) is one compo- 
nent of the deterministic limit-cycle solution [vn(f), Wo(f), ao(t)] 
that we calculated numerically. 

In Figure 3B we demonstrate that all possible correlation pat- 
terns can be realized in the GIF model and that the predicted 
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FIGURE 3 | Possible correlation patterns and corresponding PRCs (solid 
lines: theory, symbols: simulations of Equation (2)). For the adapting LIF 
model (A), ■& < 1 and only three qualitative different cases are possible. The 
adapting GIF model (B) exhibits the full repertoire of correlation patterns 
because the PRC can be partly negative and 0 can attain values from its 
entire physically meaningful interval [— 1/a, 1 /a]. The value of 0 and hence 
the type of correlation pattern is set by the integral over the weighted PRC 
Z(f) = Z(f)e _i » shown in left panels. LIF parameters: D= 0.1, 
x a = 2, (i) |i = 20, A = 10, (ii) n = 20, A = 4.47, (iii) u, = 5, A = 1 . GIF 
parameters: (i) |a = 10, p = 3, x a = 10. (ii) |i = 11 .75, f> = 3, x a = 10. (iii) 
u. = 20, f3= 1.5, x a = 10. (iv) (i = 2.12, p= 1.5, x s = 1, A = 10. (v) 
Ii = 1.5, p = 1.5, x a = 1, A = 9 D = 10~ 5 . Unless stated otherwise, y = 1, 
A = 1, x w = 1.5, D = 10~ 4 , w f = 0. 



SCCs agree quantitatively well in theory and model simulations 
(for comparison, see the SCC for the LIF in Figure 3A). To each 
distinct pattern belongs a range of ■& (Figure 3, left), determined 
by the area under the weighted PRC Z(t) = ^fe~^Z(t). The 

function Z(t) (left column in Figures 3A,B) illustrates, why an 
adapting GIF neuron can show vanishing (Figure 3Biv) or even 
purely positive ISI correlations (Figure 3Bv). In case of type II 
resetting, inhibitory input can shorten the ISI because of the neg- 
ative part in the PRC; here inhibition acts like an excitatory input. 
Consequently, a short ISI will induce a stronger inhibition (adap- 
tation) that now causes a likewise short interval and results thus 
in a positive correlation between adjacent ISIs. Also, the shorten- 
ing effect of the adaption current in the early negative phase of 
the PRC can be exactly balanced by the delaying effect of the late 
positive phase of the PRC (pseudo-renewal case, in which the area 
under Z is zero). 

3. DISCUSSION 

We have found a general relation between two experimentally 
accessible characteristics: the serial interval correlations and the 
phase response curve of a noisy neuron with spike-triggered 
adaptation. The theory predicts distinct correlation patterns like 
short-range negative and oscillatory correlations that have been 



observed in experiments (Ratnam and Nelson, 2000; Nawrot 
et al., 2007) and in simulation studies of adapting neurons 
(Chacron et al., 2000; Liu and Wang, 2001). 

Beyond negative and oscillatory correlations, we have found, 
however, that resonator neurons with spike-frequency adaptation 
can exhibit purely positive ISI correlations or a pseudo-renewal 
process with uncorrelated intervals. Adaptation currents that are 
commonly associated with negative ISI correlations (Wang, 1998; 
Chacron et al, 2001; Liu and Wang, 2001; Chacron et al, 2003; 
Benda et al., 2010; Nesse et al., 2010) can thus induce a rich 
repertoire of correlation patterns. Despite the multitude of pat- 
terns, there is a universal limit for the cumulative correlations at 
high firing rates [cf. Equation (11)], which shows that the long- 
term variability of the spike train is in this limit always reduced 
in agreement with experimental studies (Ratnam and Goense, 
2004). 

Our analytical results apply to arbitrary adaptation strength 
and time scale but require that (1) the noise is weak and white, 
(2) the deterministic dynamics shows periodic firing with equal 
ISIs (i.e., a limit-cycle exists) and (3) the adaptation current is 
purely spike-triggered with (4) a single exponential decay time. 
Regarding the weak-noise assumption, we found from numerical 
simulations quantitative agreement with our theory for values of 
the coefficient of variation (CV) up to 0.4, which is, for instance, 
typical for neurons in the sensory periphery (Ratnam and Nelson, 
2000; Neiman and Russell, 2004; Vogel et al, 2005). This holds 
even in the subthreshold regime at low CVs, where the determin- 
istic system does not follow a limit cycle. In this case, T* has to be 
replaced by the mean ISI. Moreover, we found qualitative agree- 
ment even for moderately strong noise with values of the CV up 
to 0.8, which is typical for cortical non-bursting neurons in vivo 
(e.g. Figure 3 in Softky and Koch (1993)). 

In the absence of a deterministic limit-cycle, i.e., in the 
fluctuation-driven regime at high CVs, different mathemati- 
cal approaches have to be employed, such as those based on 
a hazard-function formalism (Muller et al., 2007; Nesse et al., 
2010; Schwalger and Lindner, 2010; Farkhooi et al, 2011). 
Furthermore, for some parameter sets, we also observed repeat 
periods of the deterministic system that involved multiple ISIs 
corresponding to a periodic ISI sequence with X, = X 1 + „, where 
the smallest period is « > 2. Such cases can realize bursting (Naud 
et al., 2008), which we did not consider in the present study. 
However, we expect that these parameter regimes yield interest- 
ing correlation patterns because already in the noiseless case a 
periodic ISI sequence exhibits correlations between ISIs. 

Regarding the last two assumptions, it seems that the analytical 
derivation cannot be easily extended to the cases of adapta- 
tion currents activated by the subthreshold membrane potential 
("subthreshold adaptation" Ermentrout et al., 2001; Brette and 
Gerstner, 2005; Prescott and Sejnowski, 2008; Deemyad et al., 
2012) and multiple-time-scale adaptation (Pozzorini et al., 2013). 
Ermentrout et al. (2001) have shown that the inclusion of sub- 
threshold adaptation can lead to type II PRCs, which according 
to our theory could qualitatively change the correlation patterns. 
An adaptation dynamics depending on the subthreshold mem- 
brane potential also involves a fluctuating component because v 
is noisy. According to Schwalger et al. (2010), this stochasticity 
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could contribute positive correlations. The combined effect of 
spike-triggered, subthreshold and stochastic adaptation currents 
on the sign of the SCC is not clear. 

The important cases of the fluctuation-driven regime and 
multiple-time-scale adaptation have been recently analyzed with 
respect to the first-order spiking statistics including the stationary 
firing rate as well as the mean response to time-dependent stim- 
uli (Richardson, 2009; Naud and Gerstner, 2012). The second- 
order statistics, which describes the fluctuations of the spike 
train ("neural variability," cf. Section 4.2) and which limits the 
information transmission capabilities of neurons, is however still 
poorly understood theoretically in these cases. How adaptation 
shapes second-order statistics in the cases of multiple adapta- 
tion time scales, fluctuation-driven spiking and sub-threshold 
adaptation is an interesting topic for future investigations. 

As an outlook we sketch, how our theory could be used to con- 
strain unknown physiological parameters by measured SCCs and 
PRCs. For instance, from the mean ISI we can estimate T* = ( T) . 
Furthermore, knowing pi = — A(a, — ft) as well as the ratio 
p 2 /p! = ot-rj one can eliminate and solve for a. This allows to 
estimate the unknown adaptation time constant x a = — T*/lna 
and the amplitude of the adaptation current 



— I a 
a 



dtZ(t)e 



(17) 



Although experimental PRCs are notoriously noisy (Izhikevich, 
2005), the integral over Z(f) determining our estimate of a* is less 
error-prone. Combining our approach with advanced estimation 
methods for the PRC (Galan et al., 2005), may thus provide an 
alternative access to hidden physiological parameters using only 
spike time statistics. 

4. MATERIALS AND METHODS 

4.1. PHASE-RESPONSE CURVES OF ADAPTING IF MODELS 

We use the phase-response curve Z(t') to characterize the shift of 
the next spike following a small current pulse applied at a given 
"phase" t' € [0, T*] of an ISI. More precisely, let us assume that 
the last spike occurred at time fo = 0. Then, the next spike time 
fi of the perturbed limit cycle dynamics v = fo(v, w) + |x — a + 
eS(f - f'), v(0) = 0, w(0) = w r , fl(0) = a*, 0 < t' < T* will be 
shifted by some amount &T(t' , g) = t\ — T*. The infinitesimal 
PRC can be defined as the limit 



Z(t') 



lim 

6^0 



?>T(t\ e) 



(18) 



where the sign has been chosen such that a spike advance (8T < 
0) due to a positive stimulation (€ > 0) leads to a positive PRC. 
The definition of Z(f) by the shift of the next spike differs from 
the PRC that describes the asymptotic spike shift but is equiva- 
lent to the so-called "first-order PRC," which is often measured in 
experiments (Netoff et al., 2012). 

4. 1. 1. Adjoint equation and boundary conditions 

The PRC can be computed using the adjoint method (see e.g. 
Ermentrout and Terman (2010)). To this end, the dynamics is 



linearized about the T* -periodic limit cycle solution yo(f) 
[vo(t)> wo(f). ^o(f)]- The linearized limit-cycle dynamics y(f) 
yo(f) + &y(t) corresponding to Equation (2) is given by 



8y = A(f)8y 



(19) 



with the Jacobian matrix 



A(t) 



/ 3/0 _3/b_ 
dv dwi 

T -ii/L 

1 dv 1 dw\ 



afo 

"^1 9wn ^ 



-i \ 



To 

0 



X N dv X N dwi 



t -iMl o 



(20) 



evaluated at v = Vo(t), w = wn(f)- The linear response 
of the ISI to perturbations of the limit-cycle dynam- 
ics in an arbitrary direction is given by the vector 
Z(f) = [Z(f), Z W1 (f), . . . , Z WN (t), Z a (t)f, where the first 
component is equal to the PRC defined above. This vec- 
tor satisfies the adjoint equation Z = — A T Z (A T denotes 
the transpose of A) with the normalization condition 
v 0 (t)Z(t) + w 0 (f)Z w (f) + a 0 (t)Z a (t) = 1. The remaining 
N + 1 boundary conditions are obtained by the following con- 
sideration: On the limit cycle V, a phase <\> : V — > [0, T*] can be 
introduced in the usual way by inverting the map f \-> yo(f) and 
setting <\> = t. Because we are interested in the shift of the next 
spike, it is useful to define the isochrons (sets of equal phase) as 
the sets of all points in phase space that will lead to the same first 
spike time. Put differently, phase points belonging to the same 
isochron will have their first threshold crossing in synchrony. As a 
consequence, the threshold hyperplane defined by the condition 
v = vj is a special isochron corresponding to the phase c(> = T* . 
Note that this definition of the phase implies that the reset 
line defined by the condition v = 0, w = w r does generally not 
correspond to 4> = 0 but to positive phases if a < a* and negative 
phases if a > a*. Thus, off-limit-cycle trajectories suffer a phase 
jump upon reset. Close to the threshold, the isochrons are parallel 
to the threshold, and thus, a perturbation perpendicular to the v- 
direction does not change the phase. This insensitivity implies the 



boundary conditions Z Wl (T*) 



. Z WN (T*) = Z a (T*) = 0. 



Note that a definition of the PRC based on the asymptotic spike 
shift would require periodic boundary conditions (Ladenbauer 
etal, 2012). 

From the above considerations, it becomes clear that the PRC 
Z(t) can be computed for t s [0, T*] by solving the system 



\Zw N / 



I dv 



-13/i 
-1 3/i 



13&^ / z \ 



N dv 
--13& 



3/o T -l_3/i_ 
■> 9wjv 1 dw\ 



T -l 3/n 



7 



(21) 
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subject to the boundary conditions 



Z(T) 



1 



1 



vo(T*) f 0 (vr, w 0 (T*)) + (jl - a* + A 



Z WI (T*) = 0, 



l,...,N. 



The PRC with respect to a is determined by 



Z a 



1 



-Z a + Z{t), 



Z a (T*) = 0. 



(22) 



(23) 



(24) 



The matrix in Equation (21) is again evaluated on the limit cycle 
at v = vo(t), w = wn(f) and is therefore time-dependent. An ana- 
lytic solution of Equation (21) is possible for one-dimensional 
models with adaptation (N = 0) or general linear IF models 
although in most cases the deterministic period T* still has to be 
computed numerically. 

4. 1.2. One-dimensional case 

In the case N = 0, the PRC satisfies the equation Z = —f'(vo)Z 
with boundary condition (22). The solution is given by Equation 
(13). In order to prove Equation (14), we compute Z a (t) from 
Equation (24) yielding 



Z a (t) 



(z a (0) + J Z(t')e i df'J 



Evaluation of this expression for t = T* leads to •& = 1 + 
^-Z a (0). Finally, using the normalization condition (f(0) + [i — 

a*)Z(0) - ^Z a (0) = 1 yields Equation (14). 

4.2. RELATION BETWEEN SECOND-ORDER STATISTICS OF SPIKE 
COUNT, SPIKE TRAIN AND INTERSPIKE INTERVALS 

A stationary sequence of spike times {. . . , ti, ti+i, . . . } is 
often characterized by the statistics of the spike train x(t) = 
8(2" — f,), the spike count N(t) = J* df' x{t') or the sequence 
of ISIs { X, = f, — f,_ i } . In particular, neural variability can be 
quantified by the second-order statistics of these different descrip- 
tions as, for instance, the spike train power spectrum 

S(f) = j dx^ x (x(t)x(t + x)), (25) 



the Fano factor 



F(t) 



(N(t) 2 ) - (N(t)} 2 
(N(t)} 



(26) 



and the coefficient of variation Cy = V ((T; — {T;)) 2 )/<T,) and 
SCC pk as defined in Equation (1). These statistics are connected 
by the fundamental relationship (Cox and Lewis, 1966) (see also 
(van Vreeswijk, 2010)) 



lim F(t) 



<r ; )lim S(f) 



r 2 



k= 1 



(27) 



It shows that the summed SCC has a strong impact on 
the long-term variability of the spike train. In partic- 
ular, a negative sum yields a more regular spike train 
on long time scales than a renewal process with the 
same Cy. 
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